Models for Socio-Environmental Data
JAGS Primer Lab
December 15, 2020
Motivation for using KnitR
You can complete all the coding exercises in the JAGS Primer using simple R scripts. However, you might be wondering about how we created the JAGS Primer key you are looking at right now. We did this using Yihui Xie’s knitr package in R. This can be a highly useful tools for organizing your Bayesian analyses. Within the same document, you can:
- Describe the model and specify the joint distribution using \(\LaTeX\) and RMarkdown.
- Write the R and JAGS code for implementing your model in JAGS.
- Run the model directly in JAGS.
- Summarize the convergence diagnostics and present results.
- Add anything else pertinent to your analysis.
Best of all, knitr can produce beautiful html files as output, which can be easily shared with collaborators. We encourage you to become familiar with knitr. We recommend Karl Broman’s knitr in a nutshell as an excellent introductory tutorial.
Exercise 1: Factoring
There is no x in the posterior distribution in equation 4. What are assuming if x is absent? Draw the Bayesian network, or DAG, for this model. Use the chain rule to fully factor the joint distribution into sensible parts then simplify by assuming that \(r\), \(K\) and \(\tau\) are independent.
There is no \(x\) because we are assuming it is measured without error.
Factoring the joint using the chain rule: \[[r,K,\tau, \mathbf{y}]=[y\mid r,K,\tau][r|K,\tau][K|\tau][\tau]\]
Simplifying under the assumption of \(r,K,\tau\) are independent: \[[r,K,\tau, \mathbf{y}]=[y\mid r,K,\tau][r][K][\tau]\]
We would obtain the same result by inspecting the DAG and using the rules:
- All random variables on at the head of arrows must be on the lhs of a conditioning symbol.
- All random variables at the tails of arrows must be on the rhs of a conditioning symbol.
- All random variables at the tails of arrows with no arrows coming into them must be expressed unconditionally.
Exercise 2: Can you improve these priors?
A recurring theme in this course will be to use priors that are informative whenever possible. The gamma priors in equation 4 of the JAGS Primer include the entire number line \(>0\). Don’t we know more about population biology than that? Lets, say for now that we are modeling the population dynamics of a large mammal. How might you go about making the priors on population parameters more informative?
A great source for priors in biology are allometric scaling relationships that predict all kinds of biological quantities based on the mass organisms (Peters, 1983; Pennycuick,1992). If you know the approximate mass of the animal, you can compute broad but nonetheless informative priors on \(r\) and \(K\). This might leave out the social scientists, but I would trust the scaling of \(r\) for people if not \(K\).
In the absence of some sort of scholarly way to find priors, we can at least constrain them somewhat. There is no way that a large mammal can have an intrinsic rate of increase exceeding 1 – many values for \(r\) within gamma(.001, .001) are far large than than that and hence are complete nonsense. We know \(r\) must be positive and we can put a plausible bound on its upper limit. The only requirement for a vague prior is that its “\(\ldots\) range of uncertainty should be clearly wider that the range of reasonable values of the parameter\(\ldots\)” (Gelman and Hill, 2009, page 355), so we could use \(r\) ~ uniform(0, 2) and be sure that it would be minimally informative. Similarly, we could use experience and knowledge to put some reasonable bounds on \(K\) and even \(\sigma\), which we can use to calculate \(\tau\) as \(\tau=\frac{1}{\sigma^{2}}\).
Peters. The ecological implications of body size. Cambridge University Press, Cambridge, United Kingdom, 1983.
C. J. Pennycuick. Newton rules biology. Oxford University Press, Oxford United Kingdom, 1992.
A. Gelman and J. Hill. Data analysis using regression and multilievel / hierarchical modeling. Cambridge University Press, Cambridge, United Kingdom, 2009.
Exercise 3: Using for loops
Write a code fragment to set vague normal priors for 5 regression coefficients – dnorm(0, 10E-6) – stored in the vector b.
for(i in 1:5){
b[i] ~ dnorm(0, .000001)
}Exercise 4: Coding the logistic regression
Write R code (algorithm 3) to run the JAGS model (algorithm 3) and estimate the parameters, \(r\), \(K\) \(\sigma\), and \(\tau\). We suggest you insert the JAGS model into this R script using the sink() command as shown in algorithm 4. because this model is small. You will find this a convenient way to keep all your code in the same R script. For larger models, you will be happier using a separate file for the JAGS code. Here is the joint distribution for our logistic model again, with the priors updated from exercise 2 and \(\tau\) expressed as a derived quantity:
\[\begin{eqnarray} \mu_{i} & = & r-\frac{rx_{i}}{K}\textrm{,}\nonumber\\[1em] \tau & = & \frac{1}{\sigma^{2}}\textrm{,}\nonumber\\[1em] \left[r,K,\sigma\mid\mathbf{y}\right] & \propto & \prod_{i=1}^{n}\textrm{normal}\left(y_{i} \mid \mu_{i},\tau\right)\textrm{uniform}\left(K\mid 0,4000\right) \textrm{uniform}\left(\sigma\mid 0, 5\right) \textrm{uniform}\left(r\mid 0,2\right)\textrm{.}\nonumber\\ \end{eqnarray}\]
We use the sink() command to create a JAGS script from our joint distribution. This file is created within R and saved in the working directory. Please note that the outer set of brackets are only required when running this code within an R markdown document (as we did to make this answer key). If you are running them in a plain R script, they are not needed.
{ # Extra bracket needed only for R markdown files
sink("LogisticJAGS.R")
cat("
model {
# priors
K ~ dunif(0, 4000)
r ~ dunif (0, 2)
sigma ~ dunif(0, 50)
tau <- 1/sigma^2
# likelihood
for(i in 1:n) {
mu[i] <- r - r/K * x[i]
y[i] ~ dnorm(mu[i], tau)
}
}
",fill = TRUE)
sink()
} # Extra bracket needed only for R markdown filesThen we run the remaining commands discussed in the JAGS Primer. Note that jm is calling the JAGS script LogisticJAGS.R we just created.
#Be sure to do this sorting. It will be important for plotting later.
Logistic <- Logistic[order(Logistic$PopulationSize),]
inits = list(
list(K = 1500, r = .2, sigma = .01),
list(K = 1000, r = .15, sigma = .5),
list(K = 900, r = .3, sigma = .01))
data = list(
n = nrow(Logistic),
x = as.double(Logistic$PopulationSize),
y = as.double(Logistic$GrowthRate))
n.adapt = 1000
n.update = 10000
n.iter = 10000
jm = jags.model("LogisticJAGS.R", data = data, inits = inits, n.chains = length(inits), n.adapt = n.adapt)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 3
## Total graph size: 212
##
## Initializing modelupdate(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
MCMCsummary(zm)## mean sd 2.5% 50% 97.5% Rhat
## K 1.237502e+03 6.270432e+01 1.128602e+03 1.232806e+03 1.375907e+03 1
## r 2.006954e-01 9.686011e-03 1.813155e-01 2.006520e-01 2.196546e-01 1
## sigma 2.864351e-02 3.040275e-03 2.345873e-02 2.839627e-02 3.531073e-02 1
## tau 1.259338e+03 2.612025e+02 8.020225e+02 1.240159e+03 1.817152e+03 1
## n.eff
## K 7016
## r 7235
## sigma 15059
## tau 16330Exercise 5: Understanding coda ojects
- Look at the first six rows of the data frame.
- Find the maximum value of \(\sigma\).
- Find the mean of \(r\) for the first 1000 iterations.
- Find the mean of \(r\) after the last 1000 iterations.
- Make two publication quality plots of the marginal posterior density of
K, one as a smooth curve and the other as a histogram. - Compute the probability that
K> 1600. Hint: what type of probability distribution would you use for this computation? Investigate the the dramatically useful R functionecdf(). - Compute the probability that 1000 <
K< 1300. - Compute the .025 and .975 quantiles of
K. Hint–use the Rquantile()function.
# Make a data frame
df = as.data.frame(rbind(zm[[1]], zm[[2]], zm[[3]]))
# Look at the first six rows:
head(df)## K r sigma tau
## 1 1243.837 0.1982945 0.02899101 1189.7984
## 2 1220.880 0.1987010 0.02724243 1347.4370
## 3 1277.394 0.1869374 0.03609160 767.6931
## 4 1317.190 0.1910025 0.03256109 943.1965
## 5 1272.104 0.2131853 0.03068146 1062.3019
## 6 1267.324 0.2106143 0.03274960 932.3695# Find the mean of r for the first 1000 iterations
mean(df$r[1: 1000])## [1] 0.2009405# Find the mean of r for the last 1000 iterations
nr = length(df$r)
mean(df$r[(nr - 1000): nr])## [1] 0.2002266# Make a publication quality plot of the marginal posterior distribution of K as a smooth curve.
# More iterations would produce a smoother cruve) and as a histogram.
plot(density(df$K), main = "", xlim = c(1000, 1500), xlab = "K") 
Fig. 2. Posterior density of K.
hist(df$K, main = "", xlim = c(1000, 1500), xlab = "K", breaks = 50, freq = FALSE) 
Fig. 2. Posterior density of K.
# Find the probability that the parameter K exceeds 1600
1 - ecdf(df$K)(1600)## [1] 1e-04# Find the probability that the parameter 1000 < K < 1300
ecdf(df$K)(1300) - ecdf(df$K)(1000)## [1] 0.8510333# Compute the .025 and .975 quantiles of K
quantile(df$K, c(.025, .975))## 2.5% 97.5%
## 1128.602 1375.907Exercise 6: Summarizing coda objects
Summarize the coda output from the logistic model with 4 significant digits. Include Rhat and effective sample size diagnostics (more about these soon). Summarize the coda output for \(r\) alone.
MCMCsummary(zm, round = 4, n.eff = TRUE)## mean sd 2.5% 50% 97.5% Rhat n.eff
## K 1237.5020 62.7043 1128.6020 1232.8057 1375.9067 1 7016
## r 0.2007 0.0097 0.1813 0.2007 0.2197 1 7235
## sigma 0.0286 0.0030 0.0235 0.0284 0.0353 1 15059
## tau 1259.3378 261.2025 802.0225 1240.1589 1817.1515 1 16330MCMCsummary(zm, params = "r", round = 4, n.eff = TRUE)## mean sd 2.5% 50% 97.5% Rhat n.eff
## r 0.2007 0.0097 0.1813 0.2007 0.2197 1 7235Exercise 7: Using MCMCchains
Extract the chains for r and sigma as a data frame using MCMCchains and compute their .025 and .975 quantiles from the extracted object. Display three significant digits. Make a publication quality histogram for the chain for sigma. Indicate the .025 and .975 Bayesian equal-tailed credible interval value with vertical red lines. Overlay the .95 highest posterior density interval with vertical lines in blue. This is a “learn on your own” problem intended to allow you to rehearse the most important goal of this class: being sufficiently confident to figure things out. Load the package HDinterval, enter HDInterval at the console and follow your nose. Also see Hobbs and Hooten Figure 8.2.1.
ex = as.data.frame(MCMCchains(zm, params = c("r", "sigma")))
class(ex)## [1] "data.frame"dim(ex)## [1] 30000 2signif(apply(ex, 2, quantile, c(.025, .975)), 3)## r sigma
## 2.5% 0.181 0.0235
## 97.5% 0.220 0.0353hist(ex$sigma,xlab = expression(sigma), ylab = "Probability density", freq = FALSE, breaks = 100, main = "")
abline(v = quantile(ex$sigma, c(.025, .975)), col = "red", lwd = 2)
abline(v = hdi(ex$sigma, .95), lwd = 2, col = "blue")
Exercise 8: Plotting data and model predictions using the MCMCpstr formatting
For the logistic model:
- Plot the observations of growth rate as a function of observed population size.
- Overlay the median of the model predictions as a solid line.
- Overlay the 95% equal-tailed credible intervals as dashed lines in red.
- Overlay the 95% highest posterior density intervals as dashed lines in blue.
- What do you note about these two intervals? Will this always be the case? Why or why not?
- What do the dashed lines represent? What inferential statement can we make about relative these lines?
zm = coda.samples(jm, variable.names = c("K", "r", "sigma", "mu"), n.iter = 10000)
BCI <- MCMCpstr(zm, params = "mu", func = function(x) quantile(x, c(.025, .5, .975)))
HPDI <- MCMCpstr(zm, params = "mu", func = function(x) hdi(x, .95))
plot(data$x, data$y, pch = 19, xlab = "Population size", ylab = "Per-capita grwoth rate (1/year)" )
lines(data$x,BCI$mu[,2], typ = "l")
lines(data$x, BCI$mu[,1], lty = "dashed", col = "red")
lines(data$x, BCI$mu[,3], lty = "dashed", col = "red")
lines(data$x, HPDI$mu[,1], lty = "dashed", col = "blue")
lines(data$x, HPDI$mu[,2], lty = "dashed", col = "blue") 
Fig. 4. Median and 95% credible intervals for predicted growth rate and posterior density of K.
The intervals are exactly overlapping in this particular case. Such overlap will not always occur. Equal-tailed intervals based on quantiles will be broader than highest posterior density intervals when marginal posterior distributions are asymmetric. There is a 95% probability that the true value of the mean per-capita population growth rate falls between these lines. Not that this differs from the prediction of a new observation of population growth rate, which would have much broader credible intervals.
Exercise 9: Creating caterpillar plots with MCMCplot
Use the MCMCplot function to create caterpillar plots to visualize the posterior densities for \(r\), and \(\sigma\) using the coda object zm from earlier. Then use MCMCplot() to explore different plotting options. There are lots of these options, including ways to make the plots publication quality, show overlap with zero, etc.
MCMCplot(zm, params = "mu")
Exercise 10: Assessing convergence
Rerun the logistic model with n.adapt = 100. Then do the following:
- Keep the next 500 iterations. Assess convergence visually with
MCMCtrace()and with the Gelman-Rubin, Heidelberger and Welch, and Raftery diagnostics. - Update another 500 iterations and then keep 500 more iterations. Repeat your assessment of convergence.
- Repeat steps 1 and 2 until you feel you have reached convergence.
- Change the adapt phase to zero and repeat steps 1 – 4. What happens?
n.adapt = 100
jm.short = jags.model("LogisticJAGS.R", data = data, inits = inits, n.chains = length(inits), n.adapt = n.adapt)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 3
## Total graph size: 212
##
## Initializing modeln.iter = 500
zm.short = coda.samples(jm.short, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
par(mfrow=c(1,2))
MCMCtrace(zm.short, pdf = FALSE)
gelman.diag(zm.short)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## K 1.02 1.06
## r 1.01 1.04
## sigma 1.00 1.01
## tau 1.01 1.02
##
## Multivariate psrf
##
## 1.02heidel.diag(zm.short)## [[1]]
##
## Stationarity start p-value
## test iteration
## K failed NA 0.00968
## r passed 201 0.13492
## sigma passed 1 0.17658
## tau passed 1 0.17311
##
## Halfwidth Mean Halfwidth
## test
## K <NA> NA NA
## r passed 2.02e-01 2.09e-03
## sigma passed 2.85e-02 3.88e-04
## tau passed 1.28e+03 3.80e+01
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.558
## r passed 1 0.676
## sigma passed 1 0.547
## tau passed 1 0.587
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 1.58e+01
## r passed 2.00e-01 1.64e-03
## sigma passed 2.88e-02 2.55e-04
## tau passed 1.24e+03 2.14e+01
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.4864
## r passed 1 0.9910
## sigma passed 101 0.0854
## tau passed 101 0.0795
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 1.17e+01
## r passed 2.01e-01 1.46e-03
## sigma passed 2.85e-02 3.03e-04
## tau passed 1.27e+03 2.71e+01raftery.diag(zm.short)## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and sn.update = 500
update(jm.short, n.iter = n.update)
n.iter = 500
zm.short = coda.samples(jm.short, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
MCMCtrace(zm.short, pdf = FALSE)
gelman.diag(zm.short)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## K 1.01 1.02
## r 1.00 1.01
## sigma 1.00 1.00
## tau 1.00 1.00
##
## Multivariate psrf
##
## 1heidel.diag(zm.short)## [[1]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.627
## r passed 1 0.539
## sigma passed 1 0.393
## tau passed 1 0.357
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 2.28e+01
## r passed 2.00e-01 2.94e-03
## sigma passed 2.89e-02 4.48e-04
## tau passed 1.24e+03 3.82e+01
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.241
## r passed 1 0.395
## sigma passed 1 0.107
## tau passed 1 0.178
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.23e+03 1.57e+01
## r passed 2.01e-01 1.65e-03
## sigma passed 2.87e-02 4.11e-04
## tau passed 1.25e+03 3.49e+01
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.0824
## r passed 101 0.0941
## sigma passed 101 0.1350
## tau passed 51 0.0599
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 1.41e+01
## r passed 2.00e-01 1.72e-03
## sigma passed 2.85e-02 3.42e-04
## tau passed 1.26e+03 2.68e+01raftery.diag(zm.short)## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and sn.update = 10000
update(jm.short, n.iter = n.update)
n.iter = 5000
zm.short = coda.samples(jm.short, variable.names = c("K", "r", "sigma", "tau"), n.iter = n.iter, n.thin = 1)
MCMCtrace(zm.short, pdf = FALSE)
gelman.diag(zm.short)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## K 1 1.01
## r 1 1.00
## sigma 1 1.00
## tau 1 1.00
##
## Multivariate psrf
##
## 1heidel.diag(zm.short)## [[1]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.970
## r passed 1 0.969
## sigma passed 1 0.722
## tau passed 1 0.545
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 5.007758
## r passed 2.01e-01 0.000651
## sigma passed 2.86e-02 0.000118
## tau passed 1.26e+03 9.865859
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.1000
## r passed 1001 0.0579
## sigma passed 1 0.1397
## tau passed 1 0.0808
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 6.454694
## r passed 2.00e-01 0.000974
## sigma passed 2.86e-02 0.000103
## tau passed 1.26e+03 8.711745
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.269
## r passed 1 0.930
## sigma passed 1 0.455
## tau passed 1 0.405
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.23e+03 4.310713
## r passed 2.01e-01 0.000570
## sigma passed 2.87e-02 0.000102
## tau passed 1.25e+03 8.254917raftery.diag(zm.short)## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## K 10 12268 3746 3.27
## r 10 11096 3746 2.96
## sigma 4 5038 3746 1.34
## tau 6 6878 3746 1.84
##
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## K 8 8831 3746 2.36
## r 24 23692 3746 6.32
## sigma 5 5391 3746 1.44
## tau 6 6520 3746 1.74
##
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## K 10 12122 3746 3.24
## r 4 5299 3746 1.41
## sigma 4 4752 3746 1.27
## tau 5 5673 3746 1.51ADVANCED Exercise 11: Coding the logistic regression to run in parallel
Append R code (algorithm 5) to the script you made in exercise 4 to run the JAGS model (algorithm 2) in parallel and estimate the parameters, \(r\), \(K\) \(\sigma\), and \(\tau\). Use the proc.time() function in R to compare the time required for the sequential and parallel JAGS run. If your computer has 3 threads, try running only 2 chains in parallel when doing this exercise. If you have fewer than 3 threads, work with a classmate that has at least 3 threads.
cl <- makeCluster(3)
initsP = function() {list(K = runif(1, 10, 4000), r = runif(1, .01, 2),
sigma = runif(1, 0, 2), .RNG.name = "base::Mersenne-Twister",
.RNG.seed = runif(1, 1, 2000))}
parallel::clusterExport(cl, c("data", "initsP", "n.adapt", "n.update", "n.iter"))
ptm <- proc.time()
out <- clusterEvalQ(cl, {
library(rjags)
jm = jags.model("LogisticJAGS.R", data = data, n.chains = 1,
n.adapt = n.adapt, inits = initsP())
update(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("K", "r", "sigma", "tau"),
n.iter = n.iter, thin = 1)
return(as.mcmc(zm))
})
ParallelTime <- proc.time() - ptm
ParallelTime## user system elapsed
## 0.004 0.000 0.782stopCluster(cl)
zmP = mcmc.list(out)
MCMCsummary(zmP)## mean sd 2.5% 50% 97.5% Rhat
## K 1.237436e+03 6.282117e+01 1.128900e+03 1.232952e+03 1.374828e+03 1
## r 2.007981e-01 9.656524e-03 1.817245e-01 2.008146e-01 2.196020e-01 1
## sigma 2.861924e-02 3.008699e-03 2.349934e-02 2.832982e-02 3.515644e-02 1
## tau 1.260514e+03 2.578836e+02 8.090775e+02 1.245983e+03 1.810876e+03 1
## n.eff
## K 2374
## r 3142
## sigma 9233
## tau 9657We rerun the model sequentially and use proc.time() again for comparison.
ptm <- proc.time()
jm = jags.model("LogisticJAGS.R", data = data, inits = inits, n.chains = length(inits), n.adapt = n.adapt)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 3
## Total graph size: 212
##
## Initializing modelupdate(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("K", "r", "sigma"), n.iter = n.iter, n.thin = 1)
SequentialTime <- proc.time() - ptm
SequentialTime## user system elapsed
## 2.016 0.006 2.024Looks as if the parallel model runs 2.59 times faster. This factor should increase the more iterations you run (why?).
You might be very tempted, in the name of reproducible science, to set the seed inside each worker to the same value. What happens when you do this? Compare the results from this run to the previous run using MCMCsummary() and MCMC(trace).
cl <- makeCluster(3)
initsP = function() {list(K = runif(1, 10, 4000), r = runif(1, .01, 2),
sigma = runif(1, 0, 2), .RNG.name = "base::Mersenne-Twister",
.RNG.seed = runif(1, 1, 2000))}
parallel::clusterExport(cl, c("data", "initsP", "n.adapt", "n.update", "n.iter"))
out <- clusterEvalQ(cl, {
library(rjags)
set.seed(10)
jm = jags.model("LogisticJAGS.R", data = data, n.chains = 1,
n.adapt = n.adapt, inits = initsP())
update(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("K", "r", "sigma", "tau"),
n.iter = n.iter, thin = 1)
return(as.mcmc(zm))
})
stopCluster(cl)
zmP = mcmc.list(out)
MCMCsummary(zm)## mean sd 2.5% 50% 97.5% Rhat
## K 1.238707e+03 63.205635682 1130.7480927 1.233152e+03 1.382807e+03 1
## r 2.006284e-01 0.009562590 0.1814912 2.006711e-01 2.192591e-01 1
## sigma 2.867054e-02 0.003024429 0.0235092 2.840384e-02 3.529137e-02 1
## n.eff
## K 1713
## r 2755
## sigma 9365MCMCsummary(zmP)## mean sd 2.5% 50% 97.5% Rhat
## K 1.239659e+03 6.227385e+01 1.134518e+03 1.233489e+03 1378.0986189 NaN
## r 2.005185e-01 9.584991e-03 1.817614e-01 2.004705e-01 0.2194044 NaN
## sigma 2.859709e-02 2.896285e-03 2.358206e-02 2.843383e-02 0.0348032 NaN
## tau 1.260124e+03 2.518494e+02 8.255847e+02 1.236885e+03 1798.1932583 NaN
## n.eff
## K 3953
## r 4096
## sigma 9014
## tau 10309MCMCtrace(zm, pdf = FALSE)
MCMCtrace(zmP, pdf = FALSE)
Repeat this exercise with fixed initial values using algorithm 6.
cl <- makeCluster(3)
myWorkers <- NA
for(i in 1:3) myWorkers[i] <- word(capture.output(cl[[i]]),-1)
initsP = list(
list(K = 1500, r = .2, sigma = 1, RNG.seed = 1,
.RNG.name = "base::Mersenne-Twister"),
list(K = 1000, r = .15, sigma = .1, RNG.seed = 23,
.RNG.name = "base::Mersenne-Twister"),
list(K = 900, r = .3, sigma = .01, .RNG.seed = 255,
.RNG.name = "base::Mersenne-Twister"))
parallel::clusterExport(cl, c("myWorkers", "data", "initsP", "n.adapt", "n.update", "n.iter"))
out <- clusterEvalQ(cl, {
library(rjags)
jm = jags.model("LogisticJAGS.R", data = data, n.chains = 1,
n.adapt = n.adapt, inits = initsP[[which(myWorkers==Sys.getpid())]])
update(jm, n.iter = n.update)
zm = coda.samples(jm, variable.names = c("K", "r", "sigma", "tau"),
n.iter = n.iter, thin = 1)
return(as.mcmc(zm))
})
stopCluster(cl)
zmP = mcmc.list(out)
MCMCsummary(zmP)## mean sd 2.5% 50% 97.5% Rhat
## K 1.236654e+03 6.207641e+01 1.129370e+03 1.232326e+03 1375.7625637 1.01
## r 2.009182e-01 9.614570e-03 1.817771e-01 2.009465e-01 0.2197338 1.00
## sigma 2.863496e-02 3.043586e-03 2.351144e-02 2.840583e-02 0.0353045 1.00
## tau 1.260038e+03 2.606176e+02 8.023057e+02 1.239324e+03 1809.0121221 1.00
## n.eff
## K 1715
## r 2558
## sigma 7247
## tau 7800MCMCtrace(zmP, pdf = FALSE)
